Noise estimation in CMB time-streams and 
fast map-making. Application to the 
BOOMERanG98 data. 
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Abstract. We describe here an iterative method for jointly estimating the noise 
power spectrum from a CMB experiment's time-ordered data, together with the 
maximum-likelihood map. We test the robustness of this method on simulated 
Boomerang datasets with realistic noise. 



1 Introduction 

The map-making problem for CMB anisotropy measurements was first con- 
sidered in the context of the COBE-DMR mission 0,^3|. It was further ex- 
tended to quick algorithms for differential measurements Jl4| taking into ac- 
count 1// noise [jl5j. Given the size of the upcoming datasets (> 10 6 pixels for 
MAP) it is essential that the map-making algorithm remains quick (typically 
0(rid log Ln) where rid is the number of time-samples, and Lm the effective 
noise- filter length). Another related question, which has been pioneered by 
[^), is how to determine the noise statistical properties from the data itself. 
We present a map-making method, based on Wright's fast algorithm, to com- 
pute iteratively a minimum variance map together with an estimate of the 
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detector noise power spectrum which is needed to compute unbiased power 
spectrum estimators of the cosmological signal. 

2 Iterative map-making - Application to simulations 
2.1 Method 

We model the data stream in the following way: 

d t = PtpA p + n t (1) 

where A p is a pixelized version of the observed sky (i.e. convolved by the 
experimental beam), and nt is the detector noise after primary deconvolution 
of any filter present in the instrumental chain (e.g. bolometer time constant, 
read-out filters). We will assume here that the experiment is a total power 
measurement, i.e. that the pointing matrix Pt p contains only one non-zero 
clement per row. We now want to estimate the minimum variance map from 
this data, i.e. the map that minimizes x 2 = (d — PZ\) T N _1 (d — PA). The 
solution is given by: 

A = (P t N -1 P)~ 1 P t N _x d (2) 

A few remarks are necessary at this point. First, the matrix to be inverted is 
huge (n p i X x ripix) so that an iterative linear solver is needed. Secondly, the 
noise correlation matrix A^, 1 has to be determined from the data itself. To 
make this tractable we assume that, at least over subsets of the time-stream, 
the noise is reasonably stationary, so that the multiplication by N _1 becomes 
a convolution operator, in other words that it is diagonal in Fourier space [J 
We thus implemented the following algorithm: 



for each stationary noise subset 

• n« = d PA b) N^- 1 = (nnt)- 1 

• A<-3 +1 ') - = (P t W*P) _1 PN (j)_1 n(^ 

endfor 



To keep the algorithm as fast as possible, we took W* to be diagonal and 
constant, so that Pt\V*P is diagonal, with each element being equal to the 
number of observations per pixel, up to a multiplicative constant. The choice 

1 This is actually only approximately true since a convolution operator is a circulant 
matrix, i.e. it assumes that the time-stream has periodic boundary conditions; 
this is however a rather good approximation for a time-stream much longer than 
the effective length of the noise-filter [n0[ 
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of the numerical value of this constant is only important for the convergence 
properties of the algorithm. In this form, the algorithm is very similar to a 

Jacobi iterative solver (in which we would have (^diag{P^lSS 1 P}j instead 

of (P^W*P) 1 ). It has the additional advantage of being very easy to com- 
pute, for very similar convergence properties. Since we made the assumption 
that each noise matrix is diagonal in Fourier space, all time-domain oper- 
ations are done in Fourier space using FFTs, thus reducing the number of 
floating point operations to 0(Ln log Lm) for each subset. 

The advantages of this iterative method are obvious: it is fast {0(rid log Ln) 
operations) and cheap in memory (O(ijv) storage). However, since we try to 
estimate both the noise power spectrum and the map in a leap-frog manner, 
the convergence properties must be studied with numerical simulations. In 
particular, the stability of the algorithm is function of the number of param- 
eters that we use to describe the noise power spectrum. In other words, we 
have to assume some regularity in the latter, which translated in practice in 
fitting the amplitude of the noise in frequency band powers of varying width. 



2.2 Simulations 

To test the properties of the algorithm, we made simulations of Boomerang 
LDB time-streams with an sCDM sky, i.e. with precisely the same scanning 
strategy and a beam smoothing comparable to the 150a channel. We chose 
this channel because it was used to determine the angular power-spectrum Cg 
of the CMB anisotropics The 10 days of data are approximately divided 
in two halves corresponding to two different azimuthal scan speeds (2 degrees, 
then 1 degree per second) at a mean elevation of ~ 45 degrees. These two scan- 
speeds excited differently mechanical/thermal resonances in the instrument 
payload, and therefore were characterized by different noise power-spectra at 
low frequency; we simulated the noise time-streams for both halves according 
to the noise power spectra as measured by our method over the 150a data. 



These power spectra are shown in figure 2.2 (left panel), for both scan 
velocities. One can see in particular that for the two degree per second spec- 
trum (dotted line) the scan-synchronous harmonics are much stronger than 
in the one degree per second spectrum (where the scan-synchronous signal is 
actually dominated by the dipole emission). For a more complete description 
of these effects, as well as a precise description of the instrument, we refer 
the reader to Brendan Crill's thesis (||). 



The right panel of figure 2.2 shows the (1 degree per second) spectrum 
of the simulated (signal+noise) data (magenta) together with the recovered 
noise spectrum (black) . One can see that most of the signal power is concen- 
trated in a smooth continuum contribution between ~ 0.1 Hz and ~ 3 — 4 Hz 
(which corresponds approximately to the beam size). This plot shows that 
a naive estimation of the noise from the data time-stream would typically 
result in overestimating the noise by ~ 10 % at the degree scale. 
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Fig. 1. Left panel: Time-stream noise spectra of the two halves of the 150a channel 
data, estimated with the method described in sectionpj]. The blue line corresponds 
to one degree per second data (1 dps), the red line to 2 dps data. Note the dif- 
ference of amplitude in the scan-synchronous lines. Right panel: Spectrum of the 
(signal+noise) 1 dps simulated data (magenta) together with the recovered noise 
power spectrum (black). 



In principle, the presence of strong scan-synchronous lines in the noise 
power spectra at low frequency, as well as the overall 1/f property of the 
noise could be taken into account by the method which would down-weight 
those frequencies accordingly, but in practice, because of the phase-coherent, 
non-gaussian nature of those scan-synchronous patterns we had to high-pass 
filter the data at 0.06 x scan speed (Hz). This changes the effective pointing 
matrix P to FP where F is the high-pass filter; however, looking at Eq. 0, 
this is equivalent to leaving P unchanged, and replacing N _1 by F^N F, 
which is straightforward in Fourier space^j. 

The Boomerang data set provides a flag with each sample data, indicating 
the mode of observation (CMB scans, source scans, SZ regions, etc.) and the 
quality of the data (either good or bad, where bad data corresponds to cosmic 
ray contamination or a scan elevation change resulting in fluctuations of the 
cold plate temperature) . We thus flagged as bad in our simulations all samples 
not corresponding to good CMB scan data. In principle, the proper way to 
deal with such a data with gaps would be to construct the orthonormal noise 
eigenmodes basis with support constrained to the valid chunks, and do the 
whole analysis on this basis instead of Fourier space. However, given the 
length of the time-stream considered here (~ 5 x 10 7 time samples) it is 
completely impractical. The solution we adopted was to model those gaps as 
being samples spent on the observation of a K internal calibrator considered 
as an additional fake pixel. We had thus to make constrained realizations of 
the noise time-stream as measured by the algorithm at each iteration step in 

2 It should be noted here that any invertible filter applied to the data does lead to 
the same maximum likelihood solution, as expected. In the case of the high-pass 
filter considered here, it is not invertible, and therefore affects the noise properties 
of the resulting map by effectively projecting out some modes in the time-stream 
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those gaps, and assign their pointing matrix element outside the map to this 
fake pixel. 

Here again, the length of the time-stream comes again as a problem, since 
the constrained realizations involve in the best case (when the separation be- 
tween gaps is bigger than the noise correlation length N COT ) inverting a matrix 
of linear size N COT , and in the worse case a matrix of linear size comparable 
to the entire time-stream. Fortunately, fast algorithms have been developed 
to answer this problem (known as linear prediction). In particular, we imple- 
mented a version of Burg's algorithm (see e.g. ||) to solve this problem. 

Finally, we used HEALpix 7' pixels ||, resulting in a 138620 pixel map. 
The input simulated map (i.e. "before observation") is shown in figure 
(upper left panel), together with the corresponding sky coverage (lower left 
panel). The different "broad bands" in the coverage file correspond to distinct 
scan periods of the flight at different constant sky elevations, illustrating 
the compromise between a uniform sky coverage and a constant elevation 
scanning strategy (to avoid atmospheric gradients). The lower right panel 
shows the difference between the input map and the map recovered from the 
time-stream by the iterative method. We can see several important things in 
this error map, namely: 

• The net effect of the high-pass filtering in the time-stream is to kill the 
largest scales in the map. The fact that the pointing matrix has been 
effectively changed by the filtering has been properly taken into account 
by the algorithm (see discussion above). 

• There is no apparent residual striping in the error map, although a more 
precise way to quantify this statement may be necessary to be able to 
compare this method with alternative map-making algorithms. 

• The apparent horizontal features in the error map are related to non- 
uniform coverage (compare to coverage map in the lower- left panel), and 
are not caused by residual striping. 

Residual striping can be seen in the upper right panel of figure ^, where the 
difference of the coadded ("naive") and iterative maps is shown. We can in 
particular see that the features follow the mean direction of the scan (as 
expected from striping) and form a non-zero angle to the iso-dec direction. 

The method is not free of caveats however. As in any iterative linear 
solver, the convergence of the solution, decomposed on the (noise) kernel 
eigenmodes, is a function of the associated eigenvalue. Thus the noisiest 
pixel modes (usually at large scales since the time-stream is high-pass fil- 
tered, either by hand or as part of the algorithm if the noise is higher at low 
frequency) take (exponentially) more time to converge. The problem scales 
as the noise matrix condition number, which is a direct function of the noise 
power spectrum dynamical range and of the scanning strategy. A solution to 
this problem Q is a multi-grid version of the present algorithm, exploiting 
the hierarchical structure of the HEALPix pixelization, and rebinning the 
time-stream accordingly. W have shown that the convergence rate on large 
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scales is then drastically reduced. Another caveat is that the distribution of 
the noise power spectrum estimation is skewed, leading to a small bias in the 
noise estimation. This is similar to the "cosmic bias" discussed by [|lj in the 
case of the C/s. A possible solution would be to approximate the likelihood 
distribution of the noise spectrum estimator to correct for this bias. Finally, 
this method does not deal yet with beam asymmetry, but fast convolution 
methods pd} ] are being developed, and will be ultimately integrated, allowing 
the addition of different channels with different beams to the same map. 



3 Conclusions - Perspectives 

We described a fast iterative map-making method to simultaneously generate 
the maximum-likelihood map and the noise power spectrum from a scanning 
experiment time stream. We tested its convergence properties on realistic 
Boomerang simulations, including the effect of complex scan strategy and 
noise power spectra, resulting in a complicated noise matrix. We concluded 
that, except for the spatially largest (and ill conditioned) modes, the map 
and noise power spectra converge very quickly. This caveat can be cured 
by means of a "convergence accelerator" , where multi-grid methods provide 
a very promising solution. In addition, this method, coupled to a fast Ci 
quadratic estimator [p2p|], provide a practical way to compute the angular 
power spectrum of CMB fluctuations in the future megapixel experiments. 
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Figure 2. Simulations - Upper left: Input map - Lower left: Coverage - Upper right: Striping map (difference between iterative and 
coadded map). - Lower right: error map (iterative case). Note the different color scales used for the difference maps. 



